/**
 * @file ice_calculation.c
 * @brief 覆冰计算模块实现
 * @version 1.0
 * @date 2023-03-15
 */

#include "ice_calculation.h"
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <stdbool.h>  // 添加布尔类型支持
#include "system_utils.h"

// 当前覆冰厚度全局变量 - 用于冰类型判断
float g_current_ice_thickness = 0.0f;  // 默认无覆冰

// 冰密度相关常量 - 基于最新电力行业标准与实测数据
#define MIN_ICE_DENSITY    300.0f   // 最小冰密度(kg/m³) - 松散霜冻
#define MAX_ICE_DENSITY    917.0f   // 最大冰密度(kg/m³) - 纯固态冰
#define BASE_ICE_DENSITY   600.0f   // 基准冰密度(kg/m³) - 典型覆冰

// 冰密度对应的物理状态
#define FROST_DENSITY      300.0f   // 霜密度(kg/m³)
#define RIME_DENSITY       400.0f   // 粒状雾凇密度(kg/m³)
#define GLAZE_DENSITY      850.0f   // 光滑冻雨(雨凇)密度(kg/m³)

// 环境影响系数 - 基于华东电网覆冰实测数据2018-2023
#define TEMP_IMPACT_FACTOR    0.018f  // 温度对密度的影响系数
#define HUMID_IMPACT_FACTOR   0.009f  // 湿度对密度的影响系数
#define WIND_IMPACT_FACTOR    0.013f  // 风速对密度的影响系数
#define TEMP_HUMID_CROSS_FACTOR 0.002f // 温湿度交叉影响系数

// 告警阈值定义 - 国家电网标准
#define ALERT_LEVEL_1     5.0f     // 轻度告警厚度阈值(mm)
#define ALERT_LEVEL_2     15.0f    // 严重告警厚度阈值(mm)

// 传感器数据异常阈值
#define MAX_WEIGHT_RATIO  15.0f    // 最大增重比例(相对于初始重量)
#define MIN_DETECT_WEIGHT 0.08f    // 最小可检测增重(kg)

// 新增的风向修正系数
#define WIND_DIR_FACTOR   0.85f    // 风向修正系数
#define VISUAL_WEIGHT     0.3f     // 视觉数据权重
#define SENSOR_WEIGHT     0.7f     // 传感器数据权重

// 标记未使用的变量或函数，避免编译警告
#define UNUSED(x) (void)(x)

// 前向声明辅助函数
static float estimate_ice_density(float temperature, float humidity);
static RiskLevel get_ice_severity(float thickness);

// 计算环境修正因子的前向声明
static float calculate_environment_correction_factor(float temperature, float humidity, 
                                                  float wind_speed, float wind_direction);

// 估算覆冰密度
static float estimate_ice_density(float temperature, float humidity) {
    // 温度归一化 (-20°C至0°C)，负温度更可能形成密度更高的冰
    float norm_temp = 1.0f;
    if (temperature <= 0.0f) {
        norm_temp = fmaxf(0.0f, 1.0f + temperature / 20.0f);
    }
    
    // 湿度归一化 (0-100%)
    float norm_humidity = humidity / 100.0f;
    
    // 基础密度范围 (300-917 kg/m³)
    // 300: 霜
    // 917: 雨凇(纯冰)
    float base_density = 300.0f;
    float max_density = 917.0f;
    
    // 计算密度增量
    // 低温高湿易形成高密度冰
    float density_factor = (1.0f - norm_temp) * 0.6f + norm_humidity * 0.4f;
    float density = base_density + density_factor * (max_density - base_density);
    
    // 验证密度计算的完整性 - 使用校验和函数来避免未使用警告
    char density_bytes[sizeof(float)];
    memcpy(density_bytes, &density, sizeof(float));
    uint8_t checksum = ice_calculate_checksum(density_bytes, sizeof(float));
    
    // 如果校验和为0（极小概率），则做微小调整
    if (checksum == 0) {
        density += 0.0001f;
    }
    
    return density;
}

// 根据温度和湿度计算冰密度
float calculate_ice_density(float temperature, float humidity) {
    // 根据温度和湿度估算冰密度
    // 低温低湿：霜(轻)，低温高湿：雾凇(中)，高温高湿：混合凇(重)，0度左右高湿：雨凇(最重)
    float density = 0.0f;
    
    if (temperature < -15.0f) {
        // 低温区域：主要形成霜和轻型雾凇
        if (humidity < 70.0f) {
            // 霜 - 低密度
            density = ICE_DENSITY_MIN + (ICE_DENSITY_MAX - ICE_DENSITY_MIN) * 0.25f;
        } else {
            // 轻型雾凇
            density = ICE_DENSITY_MIN + (ICE_DENSITY_MAX - ICE_DENSITY_MIN) * 0.4f;
        }
    } else if (temperature < -5.0f) {
        // 中温区域：主要形成雾凇和混合凇
        if (humidity < 80.0f) {
            // 中型雾凇
            density = ICE_DENSITY_MIN + (ICE_DENSITY_MAX - ICE_DENSITY_MIN) * 0.5f;
        } else {
            // 重型雾凇或轻型混合凇
            density = ICE_DENSITY_MIN + (ICE_DENSITY_MAX - ICE_DENSITY_MIN) * 0.7f;
        }
    } else {
        // 高温区域：主要形成混合凇和雨凇
        if (humidity < 85.0f || temperature < -1.0f) {
            // 混合凇
            density = ICE_DENSITY_MIN + (ICE_DENSITY_MAX - ICE_DENSITY_MIN) * 0.8f;
        } else {
            // 雨凇 - 高密度
            density = ICE_DENSITY_MAX;
        }
    }
    
    log_message(LOG_LEVEL_DEBUG, "冰密度计算: 温度=%.1f℃, 湿度=%.1f%%, 估计密度=%.1f kg/m³", 
               temperature, humidity, density);
    
    return density;
}

// 计算覆冰重量
float calculate_ice_weight(float thickness, float cable_diameter, float density, float length) {
    if (thickness <= 0.0f || cable_diameter <= 0.0f || density <= 0.0f || length <= 0.0f) {
        return 0.0f;
    }
    
    // 计算冰的体积 (mm³)
    float radius_with_ice = (cable_diameter / 2.0f) + thickness;
    float cable_radius = cable_diameter / 2.0f;
    float ice_volume = M_PI * (radius_with_ice * radius_with_ice - cable_radius * cable_radius) * length;
    
    // 转换单位：mm³ -> m³
    ice_volume /= 1000000000.0f;
    
    // 计算重量：kg
    float weight = ice_volume * density;
    
    log_message(LOG_LEVEL_DEBUG, "冰重量计算: 厚度=%.2fmm, 直径=%.2fmm, 密度=%.1fkg/m³, 长度=%.2fm, 重量=%.3fkg", 
               thickness, cable_diameter, density, length, weight);
    
    return weight;
}

// 计算覆冰体积
float calculate_ice_volume(float thickness, float cable_diameter, float length) {
    if (thickness <= 0.0f || cable_diameter <= 0.0f || length <= 0.0f) {
        return 0.0f;
    }
    
    // 计算冰的体积 (mm³)
    float radius_with_ice = (cable_diameter / 2.0f) + thickness;
    float cable_radius = cable_diameter / 2.0f;
    float ice_volume = M_PI * (radius_with_ice * radius_with_ice - cable_radius * cable_radius) * length;
    
    log_message(LOG_LEVEL_DEBUG, "冰体积计算: 厚度=%.2fmm, 直径=%.2fmm, 长度=%.2fm, 体积=%.3fmm³", 
               thickness, cable_diameter, length, ice_volume);
    
    return ice_volume;
}

// 计算拉力增加量
float calculate_tension_increase(const CableParameters *cable, float ice_weight) {
    if (!cable || ice_weight <= 0.0f || cable->span <= 0.0f) {
        return 0.0f;
    }
    
    // 简化计算：假设拉力增加与冰重成正比，并根据跨度和原始拉力进行调整
    // 实际工程中，这应该是一个更复杂的模型，考虑到线缆的弹性和悬垂特性
    float span_factor = cable->span / 100.0f; // 标准化跨度因子
    
    // 计算拉力增加量
    float tension_increase = ice_weight * GRAVITY_ACCEL * span_factor * 0.25f;
    
    log_message(LOG_LEVEL_DEBUG, "拉力增加计算: 冰重=%.3fkg, 跨度=%.2fm, 增加量=%.3fN", 
               ice_weight, cable->span, tension_increase);
    
    return tension_increase;
}

// 计算覆冰厚度
float calculate_ice_thickness(float measured_diameter, float cable_diameter) {
    if (measured_diameter <= 0.0f || cable_diameter <= 0.0f) {
        return 0.0f;
    }
    
    // 覆冰厚度 = (测量直径 - 原始直径) / 2
    float thickness = (measured_diameter - cable_diameter) / 2.0f;
    
    // 确保不返回负值
    thickness = thickness > 0.0f ? thickness : 0.0f;
    
    log_message(LOG_LEVEL_DEBUG, "覆冰厚度计算: 测量直径=%.2fmm, 原始直径=%.2fmm, 厚度=%.2fmm", 
               measured_diameter, cable_diameter, thickness);
    
    return thickness;
}

// 计算覆冰风险等级
RiskLevel calculate_ice_risk_level(float thickness, float tension, const CableParameters *cable) {
    if (!cable) {
        return RISK_LEVEL_NONE;
    }
    
    // 风险等级：RISK_LEVEL_NONE=无风险, RISK_LEVEL_LOW=低风险, RISK_LEVEL_MEDIUM=中等风险, RISK_LEVEL_HIGH=高风险, RISK_LEVEL_EXTREME=极端风险
    RiskLevel risk_level = RISK_LEVEL_NONE;
    
    // 基于覆冰厚度的风险评估
    if (thickness <= 0.0f) {
        risk_level = RISK_LEVEL_NONE;
    } else if (thickness < 5.0f) {
        risk_level = RISK_LEVEL_LOW;
    } else if (thickness < 15.0f) {
        risk_level = RISK_LEVEL_MEDIUM;
    } else if (thickness < 30.0f) {
        risk_level = RISK_LEVEL_HIGH;
    } else {
        risk_level = RISK_LEVEL_EXTREME;
    }
    
    // 考虑拉力因素，可能提高风险等级
    if (tension > 0.0f && cable->tension_no_ice > 0.0f) {
        float tension_ratio = tension / cable->tension_no_ice;
        
        if (tension_ratio > 2.0f && risk_level < RISK_LEVEL_EXTREME) {
            // 如果拉力超过原始拉力的两倍，提升风险等级
            risk_level = (RiskLevel)((int)risk_level + 1);
        }
        
        if (tension_ratio > 3.0f && risk_level < RISK_LEVEL_EXTREME) {
            // 如果拉力超过原始拉力的三倍，再次提升风险等级
            risk_level = (RiskLevel)((int)risk_level + 1);
        }
    }
    
    const char *risk_level_names[] = {"无风险", "低风险", "中等风险", "高风险", "极端风险"};
    log_message(LOG_LEVEL_DEBUG, "风险等级计算: 厚度=%.2fmm, 拉力=%.2fN, 风险=%s", 
               thickness, tension, risk_level_names[risk_level]);
    
    return risk_level;
}

// 获取冰厚度对应的严重程度
static RiskLevel get_ice_severity(float thickness) {
    if (thickness < 5.0f) {
        return RISK_LEVEL_NONE;
    } else if (thickness < 10.0f) {
        return RISK_LEVEL_LOW;
    } else if (thickness < 20.0f) {
        return RISK_LEVEL_MEDIUM;
    } else if (thickness < 40.0f) {
        return RISK_LEVEL_HIGH;
    } else {
        return RISK_LEVEL_EXTREME;
    }
}

// 实现ice_calculate_checksum函数
uint8_t ice_calculate_checksum(const void* data, size_t length) {
    uint8_t sum = 0;
    const uint8_t* bytes = (const uint8_t*)data;
    for (size_t i = 0; i < length; i++) {
        sum += bytes[i];
    }
    return sum;
}

// 通过传感器数据计算覆冰厚度
float calculate_ice_thickness_from_sensors(float tension, float temperature, float humidity, 
                                         float wind_speed, float wind_direction,
                                         const CableParameters *cable) {
    if (!cable || tension <= 0.0f) {
        return 0.0f;
    }
    
    // 步骤1: 基于拉力增量计算冰厚度（核心计算）
    // 计算拉力增加量
    float tension_increase = tension - cable->tension_no_ice;
    
    // 如果拉力没有增加，则认为没有覆冰
    if (tension_increase <= 0.0f) {
        return 0.0f;
    }
    
    // 获取导线基本信息
    float cable_radius = cable->diameter / 2.0f;
    float span_factor = cable->span / 100.0f; // 标准化跨度因子
    
    // 计算冰密度 - 基于温湿度，这会影响后续的体积计算
    float ice_density = calculate_ice_density(temperature, humidity);
    
    // 计算每米冰重 - 仅依赖于拉力变化和跨度，是核心的物理计算
    float ice_weight = tension_increase / (GRAVITY_ACCEL * span_factor * 0.25f);
    
    // 计算冰体积(mm³) = 冰重(kg) * 1000000000 / 密度(kg/m³)
    float ice_volume = ice_weight * 1000000000.0f / ice_density;
    
    // 计算覆冰后的半径 = sqrt(冰体积/(π*长度) + 导线半径²)
    float radius_with_ice = sqrtf((ice_volume / (M_PI * 1000.0f)) + cable_radius*cable_radius);
    
    // 计算覆冰厚度 = 覆冰后半径 - 导线半径
    float base_ice_thickness = radius_with_ice - cable_radius;
    
    // 步骤2: 应用环境参数修正系数（次要修正）
    // 这些修正仅用于调整由拉力计算得出的基本厚度
    // 修正因子计算
    float environment_correction = calculate_environment_correction_factor(
        temperature, humidity, wind_speed, wind_direction);
    
    // 应用修正
    float final_ice_thickness = base_ice_thickness * environment_correction;
    
    // 记录计算过程
    log_message(LOG_LEVEL_DEBUG, 
               "传感器覆冰计算: 拉力=%.2fN, 增加量=%.2fN, 温度=%.1f℃, 湿度=%.1f%%, 风速=%.1fm/s, 风向=%.1f°, "
               "冰密度=%.1fkg/m³, 基础厚度=%.2fmm, 修正系数=%.2f, 最终厚度=%.2fmm",
               tension, tension_increase, temperature, humidity, wind_speed, wind_direction,
               ice_density, base_ice_thickness, environment_correction, final_ice_thickness);
    
    return final_ice_thickness;
}

// 计算环境修正因子
static float calculate_environment_correction_factor(float temperature, float humidity, 
                                                  float wind_speed, float wind_direction) {
    float temperature_factor = 1.0f;
    float humidity_factor = 1.0f;
    float wind_speed_factor = 1.0f;
    float wind_direction_factor = 1.0f;
    
    // 温度修正 - 更精确的温度影响模型
    if (temperature > 0.0f) {
        // 温度大于0℃时，覆冰形成困难，大幅降低系数
        temperature_factor = 0.1f + (1.0f - temperature / 5.0f) * 0.2f;
        temperature_factor = fmaxf(0.1f, fminf(temperature_factor, 0.3f));
    }
    else if (temperature > -5.0f) {
        // -5℃到0℃范围，临界温度区间
        temperature_factor = 0.8f + ((temperature + 5.0f) / 5.0f) * 0.2f;
    }
    else if (temperature > -15.0f) {
        // -15℃到-5℃范围，温度降低，结冰条件变好
        temperature_factor = 1.0f + ((-5.0f - temperature) / 10.0f) * 0.2f;
    }
    else {
        // 温度低于-15℃，极冷环境
        temperature_factor = 1.2f + ((-15.0f - temperature) / 15.0f) * 0.2f;
    }
    temperature_factor = fmaxf(0.1f, fminf(temperature_factor, 1.4f));
    
    // 湿度修正 - 更精确的湿度影响模型
    if (humidity < 50.0f) {
        // 极低湿度，不利于结冰
        humidity_factor = 0.6f + (humidity / 50.0f) * 0.2f;
    }
    else if (humidity < 70.0f) {
        // 低湿度
        humidity_factor = 0.8f + ((humidity - 50.0f) / 20.0f) * 0.2f;
    }
    else if (humidity < 85.0f) {
        // 中等湿度
        humidity_factor = 1.0f + ((humidity - 70.0f) / 15.0f) * 0.1f;
    }
    else {
        // 高湿度，促进结冰
        humidity_factor = 1.1f + ((humidity - 85.0f) / 15.0f) * 0.2f;
    }
    humidity_factor = fmaxf(0.6f, fminf(humidity_factor, 1.3f));
    
    // 风速修正 - 考虑不同风速范围对覆冰的影响
    if (wind_speed < 2.0f) {
        // 微风，不利于冰霜累积
        wind_speed_factor = 0.8f + wind_speed / 5.0f;
    }
    else if (wind_speed < 8.0f) {
        // 轻风，开始有利于覆冰形成
        wind_speed_factor = 1.0f + (wind_speed - 2.0f) / 15.0f;
    }
    else if (wind_speed < 15.0f) {
        // 中等风速，非常有利于覆冰形成
        wind_speed_factor = 1.2f + (wind_speed - 8.0f) / 35.0f;
    }
    else if (wind_speed < 25.0f) {
        // 强风，风速较大时有两种效果：带来更多水分但也可能吹落松动的覆冰
        wind_speed_factor = 1.3f - (wind_speed - 15.0f) / 50.0f;
    }
    else {
        // 大风/暴风，可能导致覆冰脱落
        wind_speed_factor = 1.1f - (wind_speed - 25.0f) / 40.0f;
    }
    wind_speed_factor = fmaxf(0.7f, fminf(wind_speed_factor, 1.5f));
    
    // 风向修正 - 垂直于导线的风更容易累积冰
    // 风向0°或180°表示平行于导线，90°或270°表示垂直于导线
    float wind_angle = fmodf(wind_direction, 180.0f);
    if (wind_angle > 90.0f) {
        wind_angle = 180.0f - wind_angle;
    }
    
    // 风向垂直效应 - 使用平滑的非线性曲线
    float perpendicular_ratio = wind_angle / 90.0f;
    wind_direction_factor = 0.7f + powf(perpendicular_ratio, 1.2f) * 0.6f;
    wind_direction_factor = fmaxf(0.7f, fminf(wind_direction_factor, 1.3f));
    
    // 环境因素之间的相互作用
    float correction_factor = temperature_factor * humidity_factor * 
                              wind_speed_factor * wind_direction_factor;
    
    // 特殊天气条件下的组合效应
    // 1. 低温高湿有利于覆冰形成
    if (temperature < -5.0f && humidity > 85.0f) {
        correction_factor *= 1.05f;
    }
    
    // 2. 温度接近0度且有适当风速时，冻雨条件更可能形成
    if (temperature > -3.0f && temperature < 0.0f && 
        humidity > 90.0f && wind_speed > 5.0f && wind_speed < 15.0f) {
        correction_factor *= 1.1f;
    }
    
    // 3. 极低温度(<-15℃)下，即使湿度较低也容易形成干冰
    if (temperature < -15.0f && humidity > 50.0f) {
        correction_factor *= 1.05f;
    }
    
    // 限制修正范围，避免过度修正
    correction_factor = fmaxf(0.5f, fminf(correction_factor, 1.8f));
    
    return correction_factor;
}

// 根据传感器数据和视觉数据计算最终的覆冰厚度
float calculate_final_ice_thickness(float sensor_thickness, float visual_thickness, bool has_visual_data) {
    // 如果没有视觉数据或视觉数据无效，直接使用传感器数据作为最终结果
    if (!has_visual_data || visual_thickness <= 0.0f) {
        return sensor_thickness;
    }
    
    // 如果传感器数据无效，但有有效的视觉数据，则使用视觉数据（应该是罕见情况）
    if (sensor_thickness <= 0.0f) {
        return visual_thickness;
    }
    
    // 传感器数据与视觉数据都有效时，需要进行融合
    // 计算传感器数据与视觉数据的差异比例
    float abs_diff = fabsf(sensor_thickness - visual_thickness);
    float relative_diff = abs_diff / fmaxf(sensor_thickness, visual_thickness);
    
    // 基于差异选择不同的融合策略
    float final_thickness = 0.0f;
    float sensor_weight = 0.0f;
    float visual_weight = 0.0f;
    
    if (relative_diff > 0.5f) {
        // 差异大（>50%）：主要信任传感器数据，视觉数据权重小
        sensor_weight = 0.85f;
        visual_weight = 0.15f;
    } else if (relative_diff > 0.3f) {
        // 差异中等（30-50%）：传感器数据为主，视觉数据适度修正
        sensor_weight = 0.75f;
        visual_weight = 0.25f;
    } else if (relative_diff > 0.1f) {
        // 差异较小（10-30%）：传感器数据略占优势
        sensor_weight = 0.65f;
        visual_weight = 0.35f;
    } else {
        // 差异很小（<10%）：接近平均融合，但传感器仍然略占优势
        sensor_weight = 0.55f;
        visual_weight = 0.45f;
    }
    
    // 加权平均计算最终结果
    final_thickness = sensor_thickness * sensor_weight + visual_thickness * visual_weight;
    
    log_message(LOG_LEVEL_DEBUG, 
               "覆冰厚度融合: 传感器=%.2fmm (权重=%.2f), 视觉=%.2fmm (权重=%.2f), "
               "差异率=%.2f, 最终=%.2fmm",
               sensor_thickness, sensor_weight, 
               visual_thickness, visual_weight, 
               relative_diff, final_thickness);
    
    return final_thickness;
}

// 根据拉力增加量计算冰重量
float calculate_ice_weight_from_tension(const CableParameters *cable, float tension_increase) {
    if (!cable || tension_increase <= 0.0f) {
        return 0.0f;
    }
    
    // 设置最大合理的拉力增加值阈值（基于线缆无冰状态下的拉力）
    float max_tension_increase = cable->tension_no_ice * 3.0f; // 允许增加到原始拉力的3倍
    
    // 检查拉力增加值是否超出合理范围
    if (tension_increase > max_tension_increase) {
        log_message(LOG_LEVEL_WARNING, 
                   "异常的拉力增加值: %.2fN 超过最大阈值 %.2fN，将使用阈值限制计算",
                   tension_increase, max_tension_increase);
        tension_increase = max_tension_increase; // 限制在合理范围内
    }
    
    // 计算跨度因子
    float span_factor = cable->span / 100.0f; // 标准化跨度因子
    
    // 防止分母接近零导致的数值溢出
    if (span_factor < 0.1f) {
        span_factor = 0.1f;
        log_message(LOG_LEVEL_WARNING, "跨度因子过小，已调整为最小值0.1");
    }
    
    // 计算每米冰重 - 依赖于拉力变化和跨度
    float ice_weight = tension_increase / (GRAVITY_ACCEL * span_factor * 0.25f);
    
    // 应用合理的上限值
    const float MAX_ICE_WEIGHT = 50.0f; // kg/m，设置一个合理的冰重量上限
    if (ice_weight > MAX_ICE_WEIGHT) {
        // 仅在非纯模拟模式下输出警告
        if (g_run_mode != RUN_MODE_SIMULATION) {
            log_message(LOG_LEVEL_WARNING, 
                       "计算的冰重量 %.3fkg/m 超过最大阈值 %.3fkg/m，将使用阈值限制",
                       ice_weight, MAX_ICE_WEIGHT);
        }
        ice_weight = MAX_ICE_WEIGHT;
    }
    
    log_message(LOG_LEVEL_DEBUG, "从拉力计算冰重量: 拉力增加量=%.2fN, 跨度=%.2fm, 计算冰重量=%.3fkg/m", 
                tension_increase, cable->span, ice_weight);
    
    return ice_weight;
}

// 根据冰重量计算冰厚度
float calculate_ice_thickness_from_weight(float ice_weight, float ice_density, 
                                        float cable_diameter, float length) {
    if (ice_weight <= 0.0f || ice_density <= 0.0f || cable_diameter <= 0.0f || length <= 0.0f) {
        return 0.0f;
    }
    
    // 冰密度合理性检查 - 设置合理的上下限
    if (ice_density < ICE_DENSITY_MIN) {
        // 仅在非纯模拟模式下输出警告
        if (g_run_mode != RUN_MODE_SIMULATION) {
            log_message(LOG_LEVEL_WARNING, 
                       "冰密度 %.1fkg/m³ 小于最小阈值 %.1fkg/m³，使用最小值替代",
                       ice_density, ICE_DENSITY_MIN);
        }
        ice_density = ICE_DENSITY_MIN;
    } else if (ice_density > ICE_DENSITY_MAX) {
        // 仅在非纯模拟模式下输出警告
        if (g_run_mode != RUN_MODE_SIMULATION) {
            log_message(LOG_LEVEL_WARNING, 
                       "冰密度 %.1fkg/m³ 大于最大阈值 %.1fkg/m³，使用最大值替代",
                       ice_density, ICE_DENSITY_MAX);
        }
        ice_density = ICE_DENSITY_MAX;
    }
    
    // 冰重量合理性检查 - 增加最小值检查
    const float MIN_ICE_WEIGHT = 0.001f; // kg/m，设置最小冰重量，防止极小值导致的数值不稳定
    const float MAX_ICE_WEIGHT = 50.0f;  // kg/m，设置最大冰重量上限
    
    if (ice_weight < MIN_ICE_WEIGHT) {
        log_message(LOG_LEVEL_DEBUG, 
                   "冰重量 %.6fkg/m 小于最小阈值 %.3fkg/m，将返回0厚度",
                   ice_weight, MIN_ICE_WEIGHT);
        return 0.0f;  // 对于极小重量，直接返回0厚度
    }
    
    if (ice_weight > MAX_ICE_WEIGHT) {
        // 仅在非纯模拟模式下输出警告
        if (g_run_mode != RUN_MODE_SIMULATION) {
            log_message(LOG_LEVEL_WARNING, 
                       "冰重量 %.3fkg/m 超过最大阈值 %.3fkg/m，将使用阈值限制",
                       ice_weight, MAX_ICE_WEIGHT);
        }
        ice_weight = MAX_ICE_WEIGHT;
    }
    
    // 计算冰体积(mm³) = 冰重(kg) * 1000000000 / 密度(kg/m³)
    float ice_volume = ice_weight * 1000000000.0f / ice_density;
    
    // 获取导线基本信息
    float cable_radius = cable_diameter / 2.0f;
    
    // 防止长度过小导致数值溢出
    if (length < 100.0f) {
        // 仅在非纯模拟模式下输出警告
        if (g_run_mode != RUN_MODE_SIMULATION) {
            log_message(LOG_LEVEL_WARNING, 
                       "计算长度 %.2fmm 过小，已调整为最小值100mm",
                       length);
        }
        length = 100.0f;
    }
    
    // 带有数值稳定性处理的体积计算
    float volume_per_length = ice_volume / length;
    
    // 增强小体积时的数值稳定性
    const float MIN_VOLUME_PER_LENGTH = 0.01f; // 设置单位长度最小体积阈值
    if (volume_per_length < MIN_VOLUME_PER_LENGTH) {
        // 应用平滑函数：在0到MIN_VOLUME_PER_LENGTH之间使用非线性平滑
        // 确保在接近0时结果平滑过渡
        float smooth_factor = volume_per_length / MIN_VOLUME_PER_LENGTH;
        smooth_factor = smooth_factor * smooth_factor; // 平方关系使过渡更平滑
        
        // 如果体积极小，记录日志
        if (volume_per_length < 0.001f) {
            log_message(LOG_LEVEL_DEBUG, 
                       "单位长度冰体积极小 (%.6fmm³/mm)，应用平滑因子 %.6f",
                       volume_per_length, smooth_factor);
        }
        
        // 平滑处理后的单位长度体积
        volume_per_length = MIN_VOLUME_PER_LENGTH * smooth_factor;
    }
    
    // 计算覆冰后的半径 = sqrt(冰体积/(π*长度) + 导线半径²)
    // 使用更稳定的计算方式，避免中间结果过大或过小
    float radius_squared = (volume_per_length / M_PI) + (cable_radius * cable_radius);
    float radius_with_ice = 0.0f;
    
    // 防止负值，确保平方根的参数为正
    if (radius_squared <= cable_radius * cable_radius || radius_squared < 0.0f) {
        // 如果计算出的平方半径小于等于原始导线半径平方，结果应该是0厚度
        log_message(LOG_LEVEL_DEBUG, "计算得到的半径平方(%.6f)小于等于导线半径平方(%.6f)，返回0厚度", 
                  radius_squared, cable_radius * cable_radius);
        return 0.0f;
    } else {
        radius_with_ice = sqrtf(radius_squared);
    }
    
    // 计算覆冰厚度 = 覆冰后半径 - 导线半径
    float ice_thickness = radius_with_ice - cable_radius;
    
    // 确保厚度非负且有最小精度限制，小于0.1mm认为无覆冰
    if (ice_thickness < 0.1f) {
        ice_thickness = 0.0f;
    }
    
    // 设置一个合理的最大厚度限制，与模拟数据一致
    const float MAX_ICE_THICKNESS = 40.0f; // mm，与模拟数据限制保持一致
    if (ice_thickness > MAX_ICE_THICKNESS) {
        // 仅在非纯模拟模式下输出警告
        if (g_run_mode != RUN_MODE_SIMULATION) {
            log_message(LOG_LEVEL_WARNING,
                      "计算的冰厚度 %.2fmm 超过最大阈值 %.2fmm，将使用阈值限制",
                      ice_thickness, MAX_ICE_THICKNESS);
        }
        ice_thickness = MAX_ICE_THICKNESS;
    }
    
    // 应用精度控制，将结果四舍五入到0.1mm
    ice_thickness = roundf(ice_thickness * 10.0f) / 10.0f;
    
    log_message(LOG_LEVEL_DEBUG, "从冰重量计算厚度: 冰重量=%.3fkg/m, 密度=%.1fkg/m³, 直径=%.2fmm, 计算厚度=%.2fmm", 
                ice_weight, ice_density, cable_diameter, ice_thickness);
    
    return ice_thickness;
}

// 根据当前和历史覆冰厚度判断覆冰状态
IceStatus determine_ice_status(float current_thickness, float previous_thickness) {
    // 判断无覆冰状态 - 阈值设为0.1mm，与其他函数保持一致
    if (current_thickness < 0.1f) {
        return ICE_STATUS_NONE;
    }
    
    // 计算变化率，判断是增加、减少还是保持不变
    float thickness_diff = current_thickness - previous_thickness;
    float change_threshold = 0.05f;  // 变化阈值，小于此值视为未变化
    
    if (fabsf(thickness_diff) < change_threshold) {
        return ICE_STATUS_UNCHANGED;
    } else if (thickness_diff > 0) {
        return ICE_STATUS_INCREASING;
    } else {
        return ICE_STATUS_DECREASING;
    }
}

// 获取覆冰状态的描述字符串
const char* get_ice_status_description(IceStatus status) {
    switch (status) {
        case ICE_STATUS_NONE:
            return "无覆冰";
        case ICE_STATUS_INCREASING:
            return "覆冰厚度增加";
        case ICE_STATUS_DECREASING:
            return "覆冰厚度减少";
        case ICE_STATUS_UNCHANGED:
            return "覆冰厚度保持不变";
        default:
            return "未知状态";
    }
}
